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THE THEORY OF THE GRAVITATIONAL POTENTIAL 


APPLIED TO ORBIT PREDICTION 

By James C. Kirkpatrick 
Lyndon B. Johns<m Space Center 


SUMBIARY 


Because of the significance of high-order-gravitaticmal-potential terms in 
trajectory prediction, an analysis was performed to determine the magnitude of 
position and velocity vector errors associated with the geopotential function. 

The analysis included a complete derivation of the geopotential function and its 
gradient, the transformation of Laplace's equation from Cartesian to spherical 
coordinates, and the analytic solution to Laplace's equation from the transformed 
version obtained in the classical manner of separating the variables. The solu- 
tion, therefore, oqpresses the gravitational potential at a point in space in 
terms of the magnitude of its position vector and its orientation in terms of its 
geocentric latitude and east longitude. The method devised by Pines, in which 
direction cosines are used to express the orientation of a point in space, also 
is considered in the solution. 

The effect of varying the order and the degree of the potential function 
was demonstrated by plotting trajectory integration data on position and velocity 
vector differences for a single orbit and for 32 orbits. In the resxilting curves, 
the data for lower order and lower degree potential functions are compared with 
data for an eighth-order model used as a standard. The short- and long- 
duration comparison studies were performed both including and excluding the 
effects of drag. Third-body perturbations considering the Sun and the Moon 
only were included in all studies performed. 

INTRODUCTION 


The question of the importance of high-order-gra 'itational-potential per- 
turbations of satellite orbits has been raised many times. The question is not 
without merit, for the inclusion of high-order-potential logic greatly increases 
the storage and executio*> time of any numerical integration of trajectory equa- 
tions of motion. Although the works of Pines (ref. 1) and Spencer (ref. 2) 
have made it possible to minimize the burden of storage and execution time for 
any trajectory integration considering high-order-potential terms, the doubting 
still persists as to the necessity for such terms. This indecision is unfortunate, 
since it is often important to include high-order-gravitational-potential effects 
for trajectory prediction. As a result, this report has been prepared to fully 



explain not only the potential function but also its gradient, since the gravita- 
tional force, being conservative in nature, is therefore derivable from the gra- 
dient of the potential. 

To explain the potential function, it will be necessary to derive it fully. 
Once the potential function has been derived, its gradient will be obtained. The 
derivation will be performed in the classical manner; that is. from Laplace's 
equation in spherical coordinates by the method of separating the verifies . As 
a result, the solution will require that the position of a point in space be ex- 
pressed in terms of the magnitude of its position vector and that its orientation 
be expressed in terms of its geocentric latitude and east Icmgitude. The work 
of Pines and Spencer (refs. 1 and 2) is essentially a reworking of the classical 
solution in terms of direction cosines to eliminate its isolated singularities over 
the poles and its costly trigonometric functions. 

The effect of varying the order of the potential function is illustrated in 
the figures. One series of illustrations considers the position and velocity vec- 
tor differences resulting from a trajectory integration over one single orbit; the 
resxUts of lower order potential functions are compared with an eighth-order 
model as the standard. Another series of illustrations presents the same com- 
parison study extended over 32 orbits. Both comparison studies are performed 
with and without inclusion of the effects of drag. 

The transformation of Laplace's equation from Cartesian to spherical coor- 
dinates is performed in appendix A. A concise formulation of the method of 
Pines is given in appendix B. 


SYMBOLS 


A.B 





a.b 


"ij 


®1’®2*®3’®4 



C 


arbitrary constants of power series 

determinant of identity matrix formed by excluding i-th row 
and j-th column 

function defined in equation (B3) 

constant coefficients of power series 
elements of the transformation matrix 

coefficients of the gravitational potential gradient defined in equa- 
tions (B29), (B30)> (B31), and (B32), respectively 

elements of the i>tverse transformation matrix 
arbitrary constant defined in equation (12) 
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nm’ nm 


nm 




F(a,b;c;d) 
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P 

nm 


(JA) 






QnO^^^> 


q,k 

R,<I>.A 


harmonic coefficient constants 

mass coefficient function defined in equation (B7) 

mass coefficient function defined in equation (Dll) 

hypergeometric series function defined in equation (44) 
mass coefficient function defined in equation (B12) 

gravitational parameter 

variable defined in equation (86) 
function defined in equation (52) 

positive integer counter used in equation (44) ; 0 < i < k - 1 
constant used in equation (92) and appendix B; i = 
positive integer used in equation (111); i = 1, 2, 3 
unit base vectors along Cartesian coordinate axes 
integer used in equation (111); k^ = 0, 0, 1 

integer used in equation (111); k^ = 1, -1, 0 

upper limit imposed to avoid negative factorials 

arbitrary constant deuaed in equation (15) 

positive integers denoting polynomial degree and order , 
respectively 

Legendre polynomial of the first kind, defined in equation (81) 

Legendre polynomial of the first kind, defined in equation (53) 

Legendre polynomial of the second kind, defined in equation (82) 

Legendre polynomial of the second kind, defined in equation (54) 

positive integers 
spherical coordinate functions 
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R’ = 3R/dr 


R" = a^R/dr^ 

R unit vector defined in equation (B9) 

equatorial radiua of attracting^ body 



Rm(s.t).L(8.t) 

m m 
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yjia 


x,y,z 
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a,p 

X 

M(n) 

^n 
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real and imaginary parts . respectively , of the complex variable 
(s + it)*” , where i = 

magnitude of position vector of reference point 
arbitrary constant 
direction cosine; s = x/r 
coefficient defined in equation (103) 

positive integer in text equations (51) and following 
direction cosine in appendix B; t = y/r 
gravitational potential 

complex number defined in equation (B15) 

component magnitudes of posits n vector along Cartesian coor- 
dinate axes 

arbitrary constant used in equations (23) , (84) , and (91); 
a = ±im, where i = 

constant coefficients of power series 

east longitude of reference point 

function defined in equation (61) 

direction cosine; p = z/r = sin 

variable defined in equation (B4) 

geocentric latitude of reference point 
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Operator 

7 


gradient vector operator 


ANALYSIS 


Laplace's Equation in Cartesian Coordinates 

The derivation of the geopotential function is obtained as the solution of one 
of the most important partial differential equations of mathematical physics: La- 
place's equation. Laplace's equation is defined as 


or 


= 0 


(1) 


— ? — 5 
dx^ dy^ 


+ 


aV 


= 0 


( 2 ) 


2 

where V V is defined as the Laplacian of the potential function V - V(x,y ,z): 

X, y, and z are the component magnitudes of the position vector along the Carte- 
sian coordinate axes; and V is the gradient vector operator 


V 




( 3 ) 


defined in Cartesian space . with i^, and k constituting a set of unit base vec- 
tors along the coordinate axes . The scalar product of 7*7 operating on V gives 
Laplace's equation: 


7*7V 



— ? — k 

3y^ 3z^ 


= 0 


( 4 ) 


The theory of the solution of Laplace's equation is called potential theory. 
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The Solution of Laplace's Equation in Spherical Coordinates 

Laplace's equation in spherical coordinates , expressed in the form developed 
in appendix A , is 




2 av\ 


1 8 
cos <p 


(< 


8V\ 

iCOS (p + 


1 

— 2 “ 
cos (p 


d^V 


dk 


= 0 


(5) 


where r is the magnitude of the position vector to the point in question , cp is the 
geocentric latitude , and X is the east longitude . This equation can be solved an- 
alytically by the technique of separating the variables , which is done by assuming 
a solution of the form 


V(r,<p.X) =R(r)«I>(<p)A(X) 


( 6 ) 


Therefore . 


8V _8R 
8r ~ ^ 


OA 


(7> 


8V 

7<p 


= R 


8<1> 

' 


( 8 ) 


8V 

5T 


= Ril> 


8A 

5X 


a'} 


Substituting equations (7) to (9) in equation (5) and dividing by the product of the 
spherical coordinate functions R4>A yields 


1 8 
R5? 




+ 


1 8 
4> cos 9 5^ 



1 

cos^tp 


1 8^A 


= 0 


( 10 ) 


or 


1 8 / 2 8R\ _ 1 

cos 9 


8 

5^ 


>cos 9 



+ K ^ 

cos ^9 ^ 8X^. 


( 11 ) 
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It is clear that the left-hand side of equation (Al) is a function of r only . There- 
fore . if the equality is to hold , each component on the right must be equal to some 
arbitrary constant, say C. Under these conditions, equation (11) becomes 


1 a 
R5F 



= C 


( 12 ) 


which leads to the differential equation 


r^R" + 2rR' - CR = 0 


(13) 


2 2 

where R' = 3R/ar and R” = d R/3r . Solving for the arbitrary constant C in 
equation (11) gives 


-C= .-j-- ^ (cos9®^).-V- 
* cos <p 3<P V ^ 5(p/ 

Equation (14) can be separated by multiplying both sides 
and rearranging to yield 


1 3^A 


( 11 ) 


of the equation by cos (p 


2 

-C cos • . 


_ cos <p 

4 


9 / 

S9 


3«1>\ _ 1 3^A 



(15) 


2 

Equation (15) implies that A'VA is equal to a constant -m , and this relationship 
establishes the differential equation 


A" + m^A ^ 0 


(16) 


Furthermore , equation (15) yields 


-C cos tp 


cos 4> 


3_ 

3i<p 




(17) 
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Rearranging equation (17) yields the differential equaticm 


CO., 5? 

(18) 

Equati<m (18) can be written in a more concise form by letting |x = sin 9 ; 
cos 9 = d|i/d 9 . and 3^/89 can be written as 

then. 

39 _ 39 3|i 

89 d|i 59 


39 

= C08 9 g^ 

(19) 

Substituting equation (19) in equaticm (18) 3 rields 



( 20 ) 

Equation (20) leads to the differential equation 


(1 - - 2 ^ 9 ' + 9 = 0 

( 21 ) 

Laplace's equation in spherical coordinates has now been separated into the 
following three second-order differential equations . 

(1 - p^)9" - 2p9' + ^ 9 = 0 

( 21 ) 

r^R" + 2rR' - CR = 0 

(13) 

A" + m^A = 0 

(16) 
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In these equations, C and m represent arbitrary (xmstants . Thus , the solution 
to equation (5) will be formed hrom equati<m (6) when the functicms 4>(<p) , R(r) , 
and A(X) have been determined by solving equations (21), (13), and (16), 
respectively . 

Equation (21) is known as Legendre's associated differential equation. The 
simpler form of this equation, which results when m = 0, is known as Legendre's 
differential equation. Equation (13) is of the form of Cauchy's (or Euler's) linear 
differential equation, and equation (16) is the differential equation of the undamped 
harmonic oscillator . As a result , solutions to equations (13) and (16) may be ob- 
tained by assuming solutions of the form 


R(r) = Cr** 

(22) 

A(X) = Ae“^ 

(23) 


where h is a variable defined in equation (86) , A is an arbitrary constant , and 
o is an arbitrary constant , and by substituting these functions in the eqipropriate 
differential equations . The resulting arbitrary constants are then evaluated from 
the initial and boundary conditions of the problem. Equations (13) and (16) will 
be solved after the solution of equation (21) has been completed. 

The solution of Legendre's associated differential equation will be obtained 
from the solution of Legendre's differential equation , which will be obtained first . 
It will be shown that Legendre’s differential equation has a solution in terms of two 
infinite series , each multiplied by one of the arbitrary constants of the solution . 
It will be shown further that , because of the boundary conditions of the problem 
(-1 < p < 1) , one of the arbitrary constants of the solution must be zero. 

Legendre's associated differential equation with m = 0 reduces to the form 



- 2n4>’ + C4> = 


0 


(24) 


If the assumption is made that Legendre's differential equation, equation (24) , has 
a power series solution of the form 


<D = 





(25) 
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where q and k are poaitive integers, then, when this expression is substituted 
in equation (24). the following expressi(m results. 


(l - p*) ^ (q - 2k) (q - 2k - 
k=0 


2 


+ C 


E v’-“' 


= 0 


(26) 


Bquaticm (26) reduces to the form 


(q - 2k) (q - 2k - D«2k»***’**'’^ 
k=0 


I(q - 2k)(q - 2k + 1) - = 0 


(27) 


Shifting the index k to k -t- 1 in the negative term in equation (27) and setting 
o _2 equal to zero since q is still to be determined gives 

OO 

(q - 2k) (q - 2k - DOjk****'**'’* - I (q - 2k - 2) (q - 2k - 1) - = 0 

k=0 


(28) 


Solving for a 2 j^ in equation (28) gives the recursive relation 


« - (g - 2fc - 2)(q - 2k - 1) - C . 

®2k “ (q - ik)(q - 2k - 1) °2k+2 


(29) 
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Since the assumed power series solution was defined only for positive values of k 
(i.e. , 0 < k < +••) , setting k = -1 in equation (29) causes a _2 = 0. Under these 
conditions, the following relationship results. 


q(q + 1) - C = 0 


(30) 


If the constant C is set equa, to 


C =n(n + 1) 


(31) 


where n is , of necessity , a positive integer greater than or equal to zero , then 
modifying the left-hand side of equation (30) by adding and stibtracting the product 
nq permits factoring of this equation accordingly . 


(q - n) (q + n + 1) = 0 


(32) 


Equation (32) has two possible solutions . 


q = n 

(33) 

q = -(n + 1) 

(34) 


The solution to equation (24) may new be written using A and B as the 
arbitrary constants of the solution . 


<I> = A 


00 

V' n-2k 
“2k^‘ 


+ B 




(35) 


where the constants 


°2k ^2k computed from the recursive relations 


„ - (n - 2k)(n- 2k- 1) . 

°2k+2 ' “2k 


„ _ (n + 2k + 1) (n + 2k + 2) „ 

P2k+2 2(k + l)(2n + 2k + S) *^2k 


(36) 

(’7) 
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If Oq and pQ are the arbitrary constants of the solution, then, for k = 0, 


a 


n(n - 1) 

M(in - ly 


( 38 ) 


h 


(n + l)(n + 2) 
l*2(2n + 3) 


Po 


(39) 


Wh n k = 1, equations (36) and (37) give 


“4 


(n - 2)(n - 3) 
M(in - 3^ 


“2 


_n(n - l)(n - 2)(n - 3) _ 
■ l«2-4(2n - lK2n - 3> “o 


(40) 



(n + 3)(n + 4) „ 

2*2f2n + 5) h 


_ (n + 1) (n + 2) (n + 3) (n + 4) « 


(41) 


Equatia (35) can now be written in terms of equations (40) and (41) , with 
Og and Pq equal to A and B within an arbitrary constant. 


4> = 


n(n - 1) n-2 

l«:(2n - 1) ^ 


. n(n - l)(n - 2)(n - 3) n-4 
l«2*4(2n - l)(2n - 3) ^ 



+ B 





(n + l)(n + 2) 
1»2 (2n + 3) 


1 ^ (n + l)(n + 2)(n + 3)(n + 4) 1 

"H+J ^ 1.2.4<2n + 3>(2n + 5) 


+ 


] 


(42) 


Hobson (ref. 3) points out that equation (42) can be written in the more compact form 


4> = Ah"f 




- n 
T' 


1 - 2n 
—T~' 




1 „/n + ln + 2 2n + 3 l\ 


(43) 


12 



where the function F(a.b;c;d) is the h 3 rpergeometric series defined as 


F(a.b;c;d) = 1 + ^d 


. a(a + l)b(b + 1) 
P^cTcTTI 


^2 . a(a + l)(a + 2)b(b + l)(b + 2) ^3 
^ l*2*3c(c + l)(c + 2) ^ 


+ . . . 


= 1 + 


k-1 

_ J7 <0 *»«> ^ » 


k=l 


ITT 

k! J7 + i) 

i=0 


(44) 


where the positive integer counter i has limits 0 < i < k - 1 . 

Equations (42) to (44) are not very useful for computational purposes . Equa- 
tion (35) may be written in a more useful form if and are written as 


® 2”n!n! 


2"n!n» 


Pq " (2n + 1)! 


(45) 


(46) 


Substituting equations (45) and (46) in equations (38) to (41) gives 

n(n - 1) . 2n(2n - l)(2n - 2)! 

“2 l«2(2n - 1) 2nn(n - l)(n - 2)!n(n - 1)! 

(2n ~ 2) ! (47) 

2"(n - l)(n - 2)! 


„ _ (n + l)(n + 2) . 2"n!n!(2n + 2) 

^2 l*2(2h + TT~ (2n + l)!(2n + 2) 


2"(n + D! (n 2)! 
(2n + 3){ 


(48) 
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= (n - 2) (n - 3) . (2n - 2)! 

2**(n - l)(n - 2)! 


(2n- 2)t 

2"2!(n - 2)(n - 4)! 


( 49 ) 


« .. (n + 3)(n + 4) . 2"(n + 1)! (n + 2)! 
•*4 4(2n + 5) (2n + i)! 


_ 2"(n + 2)! (n + 4)! 
2(2n + 5)! 


(50) 


Equation (42) can now be written as 


♦ Qi) =A 


1 ^ (-l)*(2n - 2t)! 

„n t! (n - t) ! (n - Ih) ! 
^ t=0 


n-2t 


+ B2" 


00 

E 


(n + k)!(n + 2k)! 1 

k! (2n + 2k + S)! ^n+2k+l 


(51) 


where t is a positive integer and S is the greatest integer equal to n/2 when n 
is an even integer and the greatest integer equal to (n - l)/2 when n is odd. 
This restriction must be imposed to avoid negative factorials. 


It is clear from equation (51) that the first series will terminate or converge 
for all values of p, including complex values. It is also immediately clear that the 
second series will diverge for all values of p such that their absolute values are 
less than unity . Since for the case of the (^potential , the value of p is bounded 
such that ~1 < p = sin 9 < 1 , the second series does not satisfy the boundary con- 
ditions and, as such, must be disregarded in the solution. Thus, the solution of 
equation (24) is written as 


I(p)=APno(»0+BQ^O^»*> 


(52) 
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where and Qj^q(m) are Legendre polynomials of the first and second kind , 

of degree n and order zero , respectively . These fianctiims are defined as 




1 V' (-l)*(2n - 2t)» ..n-2t 

Z^xHn -\)\(n - iO! ^ 

^ t=0 


(53) 


QnO<»^ = 2*^ 


PO 

E 

i»-n 


(n + k)! (n + 2k)! 

kK7n + 2k'™3>r 


1 

n+2k+l 

I* 


(54) 


Therefore , the value of the arbitrary constant B in equation (52) must be zero to 
satisfy the boundary conditions . It should be said, in passing, that the conspicu- 
ous factors 2 ” and 2^ are there for the purpose of making the value or the poly- 
nomials equal to unity when n ^ 0 and p = ±1 . 

By returning to Legendre's associated differential equation (eq. (21)) and 
replacing the arbitrary constant C by its more convenient form n(n 1) . dis- 
covered in the previous solution . the following expression is obtained . 


(l - p^)<I>»’ - 2p4>' + n(n + 1) - ^ 4> = 0 (55) 

L 1 - . 


Equation (55) is Legendre's associated differential equation of degree n and order 

2 

m. A power series solution of equation (55) is inconvenient because of the (1 - p ) 
term appearing In the denominator . A more successful approach is to consider a 
solution of the form 



m/2 

M(p) 


(56) 


Substituting equation (56) in equation (55) gives a differential equation in M(p) , 
2 

with (1 - p ) as a common factor . The functions and <h" are computed 
from equation (56) to be 


<f’(p) 





-1 


M(p) +M'(p) 


.] 


(57) 
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( 58 ) 


4>"0i) = (l - - 2)(l - n*) ^ - m(l - \i^) ^]m(h) 

- 2|im (l - ji*) ^M» (^) + M" (^)| 

Substituting equations (56) to (58) in equation (55) gives 

(l - ^*)M"(n) - 2|i(m + l)M'(n) +jmn^(m - 2)(l - ji^) ^ - m + 2mii^(l - |i^) 

+ n(n + 1) - m^(l - juiKji 


M(^) = 0 


Simplifying equation (59) yields 


(l - |i*)M"(n) - 2n(m + l)M’(|x) + [n(n + 1) - m(m + 1 )]M(h) = 0 1 
(l - |i^)M"(|i) - 2p(m + l)M’(p) + (n - m)(n + m + l)M(ji) = 0 ] 

It is possible to obtain a power series solution of equation (60) by assuming 
the solution 


M(^) = ^ 


wherein a„. is a constant coefficient. Substituting equation (61) in equation (60) 
gives 


00 oo 

(l - ^^) (q - 2k) (q - 2k - - 2ji(m + 1)^ {q- 2 k)a 2 i^p**'^^'^ 


QS 

(n - m) (n + m + 1) ®2k^^ ® 
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^ (q - 2k)(q - 2k - KQ “ 2k)(q - 2k + 2m + 1) 

k=0 k=0 


- (n - m) (n + m + 1)1 a 2 k*^** = 0 

(63) 


Shifting the index in the second summation term from k to k 1 and again requir- 
ing a _2 = 0 gives 


^ (q - 2k) (q - 2k - ^ l(q - 2k - 2) (q - 2k + 2m - 1) 

k=0 k=0 

- (n - m) (n + m + 1)1 «^2k+2^**'^^'^ " ® 

(64) 


The recursive relation may now be obtained by solving for a„j^ in equation (64) . 
Thus , 


®2k""^ 




AK. T A Ul X} yii ’ 

(q - 2k ) (q - 2k - 1) 


The numerator of equation (65) can be factored as follows . 


*2k 


= [(g - 2k - 2) 


(n - m)l [(q - 2k - 2) 
(q - 2k)lq - 2K"i) 


(n + m 1)1 


a. 


2k+2 


( 66 ) 


When k=-l, = 0 and the following solutions are possible . 


qj^ * n - m 

(67) 

q 2 = " (n + m + 1) 

(68) 
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Since for a power series solution , q must be constrained to positive integer values, 
equation (68) must be associated with a solution to be disregarded. From equation 
(67) , it is seen that the relation between n and m must always be such that 


n > m 


(69) 


The solution of equation (55) may now be written with A and B as the ar- 
bitrary constants of the solution as 


«(p) = 







-n-m-2k-li 


k=0 


(70) 


where the coefficients S 2 ^ and b 2 j^ are obtained from the recursive relations 


®2k+2 " 


(n - m - 2k) (n - m - 2k - 1) _ 
2 (R +T>(2ii - Zft - n ®2k 


(71) 


_ (n + m + 2k + 1) (n + m + 2k + 2) 
2k+2 2dc l)(2n + 2k + S) °2k 


(72) 


Equations (71) and (72) are identical to equations (36) and (37) , respectively , 
when m is equal to zero. It is interesting to note that m appears only in the nu- 
merator, bearing opposite signs in each recursive relation. Taking for a. and 
bg the following values 


(2n)! 

® 2%!(n - m)! 


(73) 



n! (n + m)! 

(!2n + i)r 


(74) 
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then for k = 0 , 


*2 = 


_ (n - m)(n - m - 1) 


TTsny 

(n - m)(n - m 
2<2n- 


1 ) 


(2n)! 


2"n!(n - m)! 


(2n - 2)! 


2"(n - l)!(n - m - 2)! 


( 75 ) 


b 


2 


_ (n + m + l)(n + m + 2)i. 

rarr55 


_ (n + m + l)(n + m + 2) 2^n! (n + m)! 

2T5H + S5 Tsmyr” 


_ 2"(n + l)!(n m + 2)! 

^5H"r3T! 


(76) 


When k = 1, 


a 


4 


(n - m - 2) (n - m - 3) _ 

ZTOiT-"35 ®2 


(n - m - 2)(n - m - 3) 
T»2(2h - 5) 


(2n - 2)! 

2"(n - D! (n - m - 2)! 


(2n - 4) ! 

2"2(n - 2)!(n - m - 4)! 


(77) 


b 


4 


_ (n + m + 3) (n + m + 4) V 

ZTOiTTS) ^2 


(n + m + 3)(n + m + 4) 2*^(n + 1)! (n + m + 2)! (2n + 4) 

z^(2!T+T) csiTTsyi TSimy 


2”(n + 2)!(n + m + 4)! 

sT^irnrr! 


(78) 
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The solution ;or ^ (^) can now be written as 


4 >(^) 




1 ^ (-l/(2n-2t)! n-m-2t 

«n f ^ t! (n - t) ! (n - m - 4t) ! 

^ t=0 


+ 




m 





(n + k)! (n + m + 2k)! 
kf(2n + 2k+ 3)! 




^n+m+ 


( 79 ) 


where S is the greatest integer equal to (n - m)/2 when (n - m) is an even in- 
teger and the greatest integer equal to (n - tn - l)/2 when (n - m) is odd. This 
restriction is imposed to avoid negative factorials . 

Equation (79) reduces to equation (SI) when m = 0. As in the case of equa- 
tion (SI) , the first series terminates or converges for all positive values of n and 
m such that n > m regardless of the value of The second series converges 
only for values ot p with a magnitude greater than \inity . As a result , the series 
must be disregarded in this solution , as the value of p is bounded such that 
-1 < p = sin 9 < -<-1 . Therefore, the solution to Legendre's associated difi'erential 
equation may be written as 




( 80 ) 


whe. 8 associated Legendre functions of the first and 

second kind , respectively , each of degree n and order m , defined as 


m 


P - (\ - V (-l)^(2n - 2t)» n 

*^nm^^^ ' „n iU^t!(n - t)!(n - m - 20! 


-m-2t 


(81) 


t=0 


m 


^ .2r «n (n + k)I(n + m + 2k). 1 

^nm^‘*^ - / 2 2^ UWTWTW. ' n+m+ST+l 


k=0 


The arbitrary constant B in equation (80) must be set equal to zero to satisfy the 
boundary conditions of the geopotential. 
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Equation (13) written with C = n(n + 1) becomes 

r^R” ! 2rR’ - n(n + 1)R = 0 
Substituting^ equation (22) in equaticm (83) gives 

oh(h - + 2ohrr*'"’ - n(n + l)or*' = 0 

The characteristic equation that results from equation (84) is 

h(h - 1) + 2h - n(n + 1) = 0 
h^ + h - n(n + 1) = 0 

By applying the quadratic formula, the values of h become 

h = -|±(n*y 


( 83 ) 


(84) 


(85) 


( 86 ) 


or 


h 


1 


= n 


(87) 


hg = -(n + 1) 


( 88 ) 


The solution of equation (83) becomes 


R = C^r" + C2r"^"'^^^ 


(89) 


Since the potential function V is expected to vanish as r , the arbitrary con- 
stant must be set equal to zero. The solution reduces to the form 


R = Cr 


-(n+1) 


(90) 
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The soluti<m to equation (16) is obtained by substituting equation (23) in 
equation (16) . This gives 


Ae 


aX 


* m*) = 


( 91 ) 


which gives, for the solution of the characteristic equation, a = ±im. where i = 
The solution of equation (16) becomes 


A(X) = Ce^ + Se“^^ 


(92) 


where C and S are arbitrary constants of the solution . By ignoring the imagi- 
nary terms . equation (92) can be simplified to 


A(X) = cos(mX) + sin(mX) 
m m 


(93) 


After substituting equations (80) , (89) , and (93) in equation (6) , the com- 
plete solution to Laplace's equation in spherical coordinates becomes 


(94) 


By setting and B equal to zero, as discussed earlier, the equation for the 
geopotential function reduces to 


V(r .q, .X) = 2^ 2^ cos(mX) + sin(mX) (95) 

n=0 m=0 ^ ^ 


The constants C_ and have been written as C and S to incorporate 

m m nm nm 

all the constants of the solution . These values are determined experimentally and 
may be found in normalized form in reference 4 and in unnormalized form in ref- 
erence 5 . 

Equation (95) has units of reciprocal length ’ ’ be written to reflect units 

of kinetic energy squared so that its gradient < ;h .. of acceleration by 
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introducing the gravitational parameter g^ » which has units of length cubed per 

time squared . and the normalizing parameter , usually the equatorial radius of 
the body. Thus. 


V(r V -£(^1 


n 


The harmonic coefficient constants C _ and S _ require additional dis- 

nm nm 

cussion . The impossibility of determining values for produced when m = 0 

is of no consequence since sin(0*X) = 0. The values of are a special set, 

known as the zonal harmonic coefficients, and are entirely dependent on latitude. 

For values of n and m greater than zero but such that n > m , these constants 
are given the name "tesseral harmonic coefficients" and are dependent on both lati- 
tude and longitude . When n = m , these coefficients are given the name "sectoral 
harmonic coefficients" and are entirely dependent on longitude . 


The Gradient of the Geopotential Function 

The gravitational force, being conservative in natwe, can be obtained from 
the gradient of the geopotential function. To obtain the gradieiit of the geopotential 
function, equation (96) can be written as 


V(r,.p.X) = ^ ^ (^) P^(sin <p)[< 

n=l m=0 


P«™(sin <p) cos(mX) + sin(mX) 


nm 


nm 


(97) 


since sin(0»X) = 0 and (R^/r)® = cos(0»X) = Pqq = 1, as may be seen from equa- 
tion (81) . The term (S£CQQ)/r is the gravitational potential associated with an 
unperturbed inverse square law of force, and the constant Cqq = ±1, depending on 
whether the force is one of attraction or repulsion, respectively. The summation 
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terms form that part oi the potential associated with the disturbing acceleration . 
From equatkm (97) , 


vv = 


where 


[’ ^00 * E (^) E »>[' 

^ I n=l m=0 


1 -2. 

r 

ir«* 

r 

m 

E 


n=l 

\ 



cos(mX) + sin(mX)l 


[©"] 


n 


2 *’nm<®^ ^ ^nm 


m=0 


^ E (["nm*'*' fj] [^mn 

m=0 ' 


+ sin(mX)| 
nm j 




+ P„„(8in 9)|V|C^ cos(mX) + sin(mX) 


ffi 


(98) 


»E 9^ ^ * yi * ^ 
~~ = T 


r r 




(99) 


(100) 


sin 9 = z/r 


( 101 ) 


7P (sin 9 ) = 
nm 


••• 

{-(s)f t’-«r * 

*■ t=0 


(102) 
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(IC3) 


(-iy(2n - 2t)! 


2"t!(n - t)!(n - m - 2t)! 




vfc _ co8(mX) + 8in(mX.)1 = mf-C 8in(mX) + S _ co8(mX)lvX 

I nm lun j i nm nm j 


X = arc t 




n > _ 3X . 3X . dX . 

5y 1+5? * 


X +y 


Substituting equation (109) in equation (108) gives 


X + y X + y 

xi- yi 

^ ~ 2~~2 
X +y 


(1C4) 

(105) 

(106) 

(107) 

(108) 

(109) 

( 110 ) 
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Therefore, by writing = x, Xg = y , and Xj = z. 




dV 


dx. 

1 



“ /R^v n n 

S (— ) 2 ‘’ran‘=ta <I»[Cot * ®nm ™(mX)] 

n=l m=0 


^ m=0 




n=l 


\ '' T /- 
y -j nmt^r, 


’(^) (?[ ■(?) ■’ ^ f * ■ * ;s) 


> ( 111 ) 


n 


Ic X« 

• [^nm ^ P^^(sin q») 


n|^Cnnj sindnX) + cos(mX)J 


/ 


where i - 1, 2, 3, kj - 0, 0, 1, and k2 = 1, -1, 0, respectively, and where C 

equals (n - m) when (n - m) is an even integer and equals (n - m - l)/2 when 
(n - m) is odd. 
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Equation (111) is useful in computing the gradient ccmiponents, but its vec- 
tor form can be written concisely as follows . 



It may be seen at once that, in the neighborhood of the poles (that is, where 
<p = ±90°) , X becomes indeterminate and the gradient of the gravitational potential 
obtained from equation (112) possesses two isolated singularities. These singular- 
ities are eliminated by Pines (ref. 1) and Spencer (ref. 2) . A condensed form of the 
work of Pines is given in appendix B . 


RESULTS AND DISCUSSION 


Both the classical method, equation (112) , and the Pines method, equation 
(B33) , have been programed and evaluated on the basis of executim time, storage 
requirements , and accuracy . The accuracy study was conducted for an integration 
time of approximately 2 days. It was found that, except in the polar rei^ons, 
both methods are identical in accuracy. However, the Pines method was approx- 
imately 15 times faster and required only one-third of the storage requirements. 
For this reason, the Pines method has been adopted as the formulatimi to use in 
any integration requiring gravitational potential perturbations, especially if both 
the gravitational potential and its gradient are required. 

A small study was conducted to show the importance of gravitational potential 
perturbations , especially those of higher orders . In this study , the position and 
velocity vector differences for various potential models of varying order and degree 
were each computed with the use of an eighth-order model as the standard . These 
differences have been plotted in figures 1 to 4. Figures 1(a) to 1(0 and 2(a) to 
2 (f> present the differences obtained over the initial orbit at 32 equidistant time 
steps in the orbit . Figures 3(a) to 3(0 and 4(a) to 4(0 present differences 
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resulting over 32 equally spaced time increments equal to the period of Keplerian 
oi^it of the initial state . The initial state of die orbit is as follows . 


Semimajor axis 

Eccentricity 

Inclination 

Longitude of the node 
Argument of periapsis 
Eccentric anomaly 
Julian date 
Keplerian period 


6629.656565 kilometers 

0.01 

0 . 7854 radian 
0.7854 radian 
0.7854 radian 
0 . 7854 radian 
2442332.17 
5372.137432 seconds 


The comparison was conducted both by considering potential perturbations <mly and 
by considering potential perturt>ations together with drag and third-body perturba- 
txon from the Sun and the Moon. The atmospheric model in reference 6 was used to 
compute density in the drag formulation, and NASA Jet Propulsion Laboratory ephem- 
eris tape data were used for the third-body and atmospheric-model perturbations . 

The legend for figures 1 to 4 is as follows . Each figure page includes six 
cvirves, which are designated with ntimbers ranging from 2 to 7. Curves marked 
with a **2’* are curves in which only second-order or second-degree terms were 
considered. Curves marked with a *’3” are curves in which both second- and thi.d- 
order or second- and third-degree terms were included. In short, the higher the 
number used to designate the curve , the more complete was the extent of the poten- 
tial function used. As a result, seventh-order curves show the lowest errors be- 
cause they approximate the eighth-order reference more closely . 

The results of the study are as foUows . High-order-gravitational-potential 
models are essential to the accuracy of a trajectory integration . This statement is 
true whether the integration is to be performed over a relatively short period , such 
as over one orbit (fig. 1(b)), or (and especially so) over longer periods (fig. 

3(b)) . These curves verify the fact that nothing is gained in accuracy by increas- 
ing the degree of the potential model without adding the corresponding tesseral 
terms . These facts are also substantiated by the velocity error curves . In fact , the 
position error curves , as expected , are virtually identical to the velocity error 
curves . Only the scale along the ordinate is different . As a result , this discussion 
of results will be confined to the position error curves because the same conclusions 
apply for the velocity error curves . 

As was stated earlier , comparisons were made with and without the effects of 
drag. In both instances, the shapes of the curves were identical except that the 
magnitude of the errors for both long- and short- duration runs was slightly higher 
when drag effects were not considered . Since the retarding force of the drag was 
absent and the effects of the Sun and the Moon are negligible, this result was ex- 
pected . The curves show additional interesting effects . 

It may be said , intuitively , that the higher the order of the model , the closer 
the agreement between the results and the higher order reference would be . The 
results of this study show this theory to be true only for the very lowest and the 
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highest order comparison curves . For the one-ort>it case , although the second- 
order model did prove to be the worst of all models compared, the fourth-ord ^r 
model proved worse than the third-order model and the fifth- . sixth- . and seventh- 
order models varied continuously. In the long-duration case, the second- and 
seventh-order models were clearly the worst and best, respectively, with the 
seventh- order model maintaining almost a constant difference with the reference 
model much unlike the rest uf the comparison curves , whi(^ showed a definite secu- 
lar trend . However , the fourth- to sixth-order curves were fairly well grouped , 
with the sixth-order curve showing a clearer and more distinct pattern . 

When the tesseral terms were removed (figs . 1(b) and 3(b)) , all curves of 
fourth degree and higher showed the same amount of error. However, when drag 
was removed, the short-duration riin (one-orbit case) seemed little affected as com- 
pared with the case when drag was present. In the long-duration run, all curves, 
regardless of the order , showed the same error , which amounted to approximately 
576 kilometers after 32 orbits . 

No comparison runs were made with the tesseral terms removed from the ref- 
erence model . However , from the results concerning the error produced when the 
tesseral terms were removed , it may be concluded that beyond the fourth degree , 
there is little to be gained in adding additional zonal terms without the benefit of 
their corresponding tesseral terms . 


CONCLUDING REMARKS 


It may be concluded from the results of this t dy that , for any trajectory in- 
tegration that extends for any appreciable length of time (especially for more than 
one orbit) , the inclusion of high-order-gravitational-potential models in the equa- 
tions of motion is an absolute requirement if any degree of accuracy is to be achieved . 
Errors of 25 kilometers and more are not uncommon if the gravitational potential 
terms are excluded . 


Lyndon B . Johnson Space Center 

National Aeronautics and Space Administration 
Houston , Texas , July 31 , 1976 
986-16-00-00-72 
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Fractions of one revolution 


(a) Eighth-order model reference with drag; lower order models with drag. 

Figure 1 Position vector differences obtained over the initial orbit. The curves 
include the following order or degree terms: (2) second only, (3) second and 
third, (4) second to fourth, (5) second to fifth, (6) second to sixth, and (7) second 
to seventh . 
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(b) Eighth-order model reference with drag; lower degree models with drag. 

Figure 1.- Continued. 




Fractions of one revolution 

(c) Eighth-order model reference with drag; lower order models without drag. 

Figure 1.- Continued. 
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Position vector error, km 



(d) Eighth-order model reference with drag: lower degree models without drag. 

Figure 1.- Continued. 
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Curve 5 
Curve 6 - P 
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Eighth-order model reference without drag; lower order models without drag. 

Figure 1.- Continued. 
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(a) Eighth- arder model reference with drag; lower order models with drag. 


Figure 2 .- Velocity vector differences obtained over the initial orbit. The curves 
include the following order or degree terms: (2) second only, (3) second and 
third, (4) seccmd to fourth, (5) second to fifth, (6) second to sixth, and (7) second 
to seventh . 
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(b) Eighth-order model reference with drag; lower degree models with drag. 

Figure 2.- Continued. 
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(c) Eighth-order model reference with drag; lower order models without drag. 

Figure 2.- Continued. 
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(d) Eighth-order model reference with drag; lower degree models without drag. 

Figure 2.- Continued. 
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(e) Bighth-ordor nuNlol roCeronoe without drag; lower «rder modela without drag. 

Figure 2.- Continued. 
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(f) Eighth-order model reference without drag; lower degree models without drag. 

Figure 2.- Concluded. 
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(a) Eighth-order model reference with drag; lower order models with drag. 

Figure 3 . - Position vector differences obtained over 32 revolutions . The curves 
include the following order or degree terms: (2) second only, (3) second and 
third, (4) second to fourth, (5) second to fifth, (6) second to sixth, and (7) second 
to seventh . 
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(b) Eighth-order model reference with drag; lower degree models with drag. 

Figure 3.- Continued. 
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(e) Eighth-order model reference without drag; lower order models without drag. 

Figure 3.- Continued. 
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(f) Eighth-order model reference without drag; lower degree models without drag. 

Figure 3.- Concluded. 
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(a) Eighth-order uiodel reference with drag; lower order models with drag. 

Figure 4.- Velocity vector differences obtained over 32 revolutions. The curves 
include the following order cr degree terms: (2) second only , (3) second and 
third, (4) second to fourth, (5) second to fifth, (6) second to sixth, and (7) second 
to seventh . 
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(b) Eighth-order model reference with drag; lower degree models with drag. 

Figure 4.- Continued. 
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(c) Eighth-order model reference with drag: lower order models without drag. 

Figure 4.- Continued. 
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(d) Eighth-order model reference with drag; lower degree models without drag. 

Figure 4,- Continued. 
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(e) Eighth-order model reference without drag; lower order models without drag 

Figure 4.- Continued. 
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(f ) Eighth-order model reference without drag; lower deg^e models without drag . 

Figure 4.- Concluded. 


APPENDIX A 


THE TRANSFORMATION OF LAPLACE'S EQUATION 
FROM CARTESIAN TO SPHERICAL COORDINATES 

To accomplish the transformation of Laplace's equati<m from Cartesian to 
spherical coordinates , it is first necessary to establish certain differential relations 
that v^xist between Cartesian and spherical coordinates . These relationships are 
established as follows . 

X = r cos <p cos X 

y = r cos 9 sin X (Al) 

z = r sin 9 

where r is the magnitude of the position vector to the point in question , 9 is the 
geocentric latitude , X is the east longitude . and x, y, and z are the component 
magnitudes of the position vector al<mg the Cartesian coordinate axes . Taking dif- 
ferentials of equation (Al) gives 

dx = cos 9 cos X dr - r sin 9 cos X d9 - r cos 9 sin X dX 

dy = cos 9 sin X dr - r sin 9 sin X d9 + r cos 9 cos X dX (A 2 a) 

dz = sin 9 dr + r cos 9 d9 

Writing equation (A 2 a) in matrix form gives 

cos 9 cos X -r sin 9 cos X -r cos 9 sin X dr 

cos 9 sin X -r sin 9 sin X r cos 9 cos X d9 (A 2 b) 

sin 9 r cos 9 0 dX 
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or, by solving for dr. dtp. and dX. 


dr 


cos 9 cos X 

-r sin 9 cos X 

-r cos 9 sin X 

-1 

dx 

d 9 

= 

cos 9 sin X 

-r sin 9 sin X 

r cos 9 cos X 


dy 

dX 


sin 9 

r cos 9 

0 


dz 


To compute the inverse of the preceding matrix . its determinant is first computed by 
expanding by minors along the third row 


/ 2 2 2 2 \ 
|det| = sin <p ^r cos X sin <p cos q> - r sin X sin <p cos qy 


- r cos 


<p^r co8^<p 


cos^X + r^cos^X sin^X^ 


2 2 - 22 
= -r sin 9 cos 9 - r cos 9 cos 9 

= -r^cos 9 


(A4) 


2 2 

making use of the trigonometric identity sin 9 + cos 9 = 1. The same identity ap- 
plies in the case of X. The inverse may now be obtained by writing the transpose 
of the matrix of co factors (i.e. , the adjoint of the matrix) and dividing each term by 
the determinant of the matrix . The cofactor of the element a., is given by 
i+i 

(-1) I A- j I , where | A^j | is the determinant of the matrix formed by excluding the 
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i-th row and j-th column of the original matrix. If b.. represents the elements of 
the inverse matrix . then ^ 


hl = 


^2 = 


^3 = 


^1 = 


*»22 = 


'* 23 " 


'>31 = 


'>32 = 


'>33 = 


^r^cos^cp cos ^^^-r^cos q>^ 
cos <p cos X 

- ^^cos^q> sin j ^-r^cos q>^ 
cos <p sin X 

^r^cos^X sin q> cos X - r^sin^X sin <p cos j^r^coa «p^ 
sin q> 

-(-r sin 9 cos 9 cos X)^^r^cos 9^ 

-(l/r)sin 9 cos X 
(r sin 9 cos 9 sin X)/j^-r^cos 9^ 

-(l/r)sin 9 sin X 

- ^ cos ^9 cos^X + r cos^9 sin^>^ j ^-r^cos 9^ 

(l/r)cos 9 

^ cos^9 sin X + r sin ^9 sin ^ ^^-r^cos 9^ 

-[(sin X)/(r cos 9)] 

- ^ cos ^9 cos X + r sin^9 cos ^l^-r^cos 9^ 

(cos X)/(r cos 9) 

(-r sin 9 cos 9 sin X cos X r sin 9 cos 9 sin X cos X) 

= 0 


(ASa) 
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Therefore, with the elements <rf the inverae matrix thus est^lished, it is possible 
to write equation (A 3 ) as follows . 


dr 


cos (p cos X cos 9 sin X sin 9 


dx 

dtp 

= 

1 . ^ 1 . ,1 

- - sin 9 cos X - - sin 9 sin X - cos 9 

r r P ^ 


dy 

dX 


sin X cos X g 

r cos 9 r cos 9 


dz 


(A 5 b) 


The required derivative relationships can now be obtained from equation (ASb) . 
These derivatives will be presented with use of the following shorthand notations . 


- ?!£.. = 

9 x ’ ^xr 3 x dr 


(A6) 


rx = cos 9 cos X 

ry = cos 9 sin X 

r„ = sin 9 
z ^ 

9 = - 4 sin 9 cos X 
^x r ^ 

9 y = - 4 ^ 

<P2 = f COS9 

_ sin X 

^ _ cos X 

X =0 

X r cos 9 

y r cos 9 

z 


r = -sin 9 cos X 

X 9 ^ 

r^j^ = -cos 9 sin X 

'yr ■ ® 

^y«p ^ ^ 

ryj^ = cos 9 cos X 

*-zr = ® 

^Z 9 = *P 

^X = ® 

1 

9 -7 sin 9 cos X 

1 

<Px<p = ■ F«>® <P cos X 

1 . 

= - sin 9 sin X 

1 . 

9 = -^ sin 9 sin X 

^ r 

*Py9 ^ " r <P sin X 

9y^ = - i sin 9 cos 

P 

9 = - — sin 9 

^Z9 r ^ 


^ _ sin X 

^ _ sin 9 sin X 

^ _ cos X 

xr 2 

r cos 9 


xX r cos 9 

^ _ cos X 

^ _ sin 9 cos X 

^ _ sin X 

yr 2 

r cos 9 

y 9 2 

r cos 9 

yX r cos 9 

X =0 

X =0 

X . = 0 

zr 

Z 9 

zX 


> (A 7 ) 
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Using the relationships given in equation (A 7 ), Laplace's equation can be trans- 
formed from Cartesian to spherical coordinates. 

From equation (Al). it is seen that x, y, and z are functions of r, 9, and 
X. Therefore, the potential function V = V(x,y ,z) is expressed in terms of r , 9 , 
and X. 


V = V r +V 9 + V.X 
X r X 9^x X X 


V =V r +V 9 +V .X 

XX xr X x9^x xX x 

V =V r +V 9 +V ,X 

yy yr y y9^y y^ y 


V =V r +V 9 +V .X 
zz zr z z9^z zX z 


(A8a) 


(A8b) 


Expanding equation (A8b) gives 

V =/Vr+Vr +V 9+V9 +V. X+V,X\r 

XX ^ rr X r xr 9r^x 9^xr Xr x X xrJ x 

* * ^r**x9 * ^99^x * ^x 

* (''rX'x * VxX * %)l'l>x ^ ''.p'f x^ * Wx * ’^x’^xx) \ <A9> 

V -/'Vr+Vr +V 9+V9 +V, X+V,x\r 

yy \ rr y r yr 9r^y 9^yr Xr y X yr^ y 

^ (^X9*’y ^ ^r*^y9 ^99^y * ^9*^F9 * * ^^^y<p) ^y 

* f rx'y " Vyx * %X^ * W ^ \xV * Vyx) \ «>»> 

V = r +V r +V 9 +V 9 + V. X +V,X \r 

zz y rr z r zr 9r^z 9^zr Xr z X zr) z 

+ ^V r+Vr^ +V 9+V9 +V. X +V.X W 

\ T9 z r *Z9 99^z 9^Z9 X9 z X z<p)^z 

+ /V .r +V r . +V .9 +V 9 . +V..X +V^X X (All) 

y rX z r zX 9X^z tp^z\ XX z X zXy z 
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Collecting terms from equations (A9) to (All) gives 

rrvx y z ) wV * y * / * 7 * / 

+ 2V (r <p +r<p +r<p\+2V.^X +rX +rX^ 
r<p^ x^x y^y z^zy rX^ x x y y z zj 

+ 2V ,/<pX +<pX + q> \ \ + y (r r +r <P_+r.X 
(pXyx X ^y y ^z z/ r\xr x X 9 ^x xX x 

+ r r +r q> +r.X +r r +r w 

yp y y<p^y yX y zr z z<p z zX z/ 

+ V /<p r + <p <P + <p ,X + o) r + cp <P + q> .X 
cp^^xr X ^x<p^x ^xX X ^yr y ^y<p y ^yX y 

+ <p r +<p <p +<p.X^+V^(X r +X q> +X^X 
^zr z ^zcp^z ^zX z) X\ xr X x<p^x xX x 

+ X r +X <p +X.X +X r +X cp +X,X^ 
yr y ytp^y yX y zr z z<p^z zX zy 


(A 12) 


where V is the gradient vector operator. Substituting the derivatives obtained 
from equation (A7) in equation (A12) gives 

2 2 2 2 2 2 
r^ + r^ + = (cos <p cos X) + (cos cp sin X) + (sin <p) 


= 1 


(A13) 


q>x^ + cpy^ + cp^^ = p sin cp cos xj ■*■ (" p sin cp sin xj + cos cpj 


(A14) 


X^+X^+X^ , 

X y z V r cos 


= (- ainX ^ ^ ^^^2 

\ r cos cp y ^r cos cp y 


1 

I T~ 

r^cbs^cp 


(A15) 


^x^x ^ ^y***y ^ *^z^z ” ^ F ^ ^ (” F ^ ^ 

<p(lcos cp) 


+ sin 
= 0 


(A 16) 
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(A17) 


+ sin (p (0) 
= 0 


<pX +q>X +<pX = (- — sin <p cos X^ (- — 
^xx^yy^zz ^ / \ ^ cos <p/ 

t (- isin <p sin x)(-52ii-) * 52EJE..0 
V r ^ /\r cos <p/ r 


= 0 


(A18) 


+ ^xk\ ~ ^ ^ *** (“ F ^) 

+ (-cos <p sin X) ^ 

^ V P cos q>/ 

= i ^sin^<p cos\ + sin^X^ (A19) 


ryj.Py + ’fy<pVy + ^y\\ ~ <p sin X) + (-sin <p sin X) | sin cp sin X^ 

+ (cos <p cos X) /-SHS-LA 
^ \r cos <p ) 

= p ^sin^<p sin\ + cos\^ (A20) 


^ ^zcp^z ^’zX^z ^ <p) + (cos cos cpj + (0) (0) 


= -cos <p 
r ^ 


(A21) 


r r 
xr X 


+ r (p + r . X 
X9^x xX X 


+ r r +r <p +r .X +r r +r tp 
yr y y<p^y yk y zr z zcp^z 


+ r . X 
zX z 


2 

F 


(A22) 
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<p r + 

Tyr y 


*Pxx\ ^ <P cos (cos <p cos X) 

+ ^- i cos <p cos r ^ ^ 

+ sin tp sin X^ (- 
I r ^ / V ^ ^ / 

‘Py<p^y ^yx\ ^ ^“T ^ ^ ^°°® <p sin X) 

+ ^- p cos <p sLi ^ ^ r ^ 

- ^-Isinpcosxyj^^'j 


(A23) 


(A24) 


"^zr^ ^ ’^z<ffz * fzk\ = ^ ^ <P) + I sin cos <p^ * (0) (0) 


(A25) 


<p r + <p <p + (p ,X + <p r + <p <p + (p .X + <p r 
^xr X ^xtp^x ^xX x ^yr y ^y9^y yX y ^zr z 


+ <p <p +Q>.X = 
^ztp^z ^zX z 


r cos (p 
(A26) 


X r +X <|. tx.x = (cos 9 cos X) * (- f 

xr X x<p'*’x xX X rcosV 


COS X \ / sin X 


1 . 0 
- sin <p cos X 


r cos <p/\ r cos 9 


sin X COS X cos 9 ^ sin X cos X sin 

2 2 

r COS 9 r cos 9 


sin X cos X 
r^co8^9 


(A27) 
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sin X) 


r cos 
^sin 


— I (cos cp si 

«P/ 


2 

sin X cos X cos <p sin X cos X sin q> sin X cos X 

_ 2 _ _ 

r cos <p r cos (p r cos (p 


(A28) 


cos (pj 


\r^'z * \<p^z * ^zX^z ^ <0){ .cos <p | + (0)(0) 

= 0 (A29) 

X r +X <p +X.X +X r +X cp +X..X +X t +X q> +X.X=0 

r't» V vm^Y YA. Y vr V VA. v zr z zfo^z za z 


ir X x<p^x xX X yr y y<p y yX y zr z z(p^z zX z 
Adding all resulting terms gives 

v^v = v„(i) . * 2%<») * ^V(O) 


.30) 



5^ 


sin (p 
- I ■ 

\ r cos <p; 


= 0 


(A31) 
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2 

If equation (A31) is multiplied through by r , there results, after some simplifica- 
tion, 


„2„_ a /_2av\ . 1 a / av\ . i aV 

= 0 ( 5 ) 


The solution of equation (5) is discussed in the section entitled "Analysis . ” 
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APPENDIX B 


THE PINES FOR»«ULATION FOR THE 
GEOPOTENTIAL FUNCTION AND ITS GRADIENT 


The classical potential function (eq. (96)) and its derivatives (eq. (Ill) or 
(112)) require that the position of a point in space be expressed in terms of geocen- 
tric latitude <p , east longitude X , and magnitude of the position vector r deter- 
mined from the component magnitudes along the Cartesian coordinate axes x , y , 
and z . Because these orientation angles require for their determination the pro- 
jection of the position vector on the equatorial reference pla; 2 . the longitude be- 
comes undeHned over the poles . What is theoretically two isolated singularities is 
actually , for computational purposes . an infinite number of singularities at each of 
two isolated regions in the immediate neighborhood of the poles , as the computer 
cannot distinguish from zero a computed noninteger number having a magnitude 
less than some specified tolerance characteristic of the computer in question . Al- 
though a polar orbit , or even a nearly polar orbit , is almost impossible to achieve 
because of a number of perturbing influences , the possibility of a singularity other 
than at the origin is eliminattid in this formulation by expressing the orientation of 
a point in space in terms of its direction cosines instead of its latitude and longitude . 
The direction cosines arc always clearly defined for any direction in space , and the 
elimination of all trigonoraetr?c functions from tiie formulaticm contributes greatly to 
decreasing execution time and storage requirements . 

The class!' il expression for the gravitational potential V exerted at a point 
in space located at a distance r from the center of the attracting body having ra- 
dius and gravitational parameter g^ is given in equation (96) as 


g /K.\ 

V(r.9,X) = j 2^ P„„{sin <p) oos(mX) * sindiij,)] (96) 

n=0 m=0 


where C ^ and S are the harmonic coefficients of the potential function and 
nm nm 

P (sin (p) represents the associated Legendre functions of the first kind, of de- 
gree n and order m. Since sin <p = z/r ~ p, where p is a direction cosine, the 
associated Legendre functions may be expressed by using Rodrigues' Theorem , 
presented here without proof. 


P 

nm 


(p) = 


m 

_ co s <p 


2"n! 


^n+m 



(Bl) 


64 



Equation (Bl) is equivalent to equation (81) . Equation (Bl) is more useftil for ana- 
lytical purposes, whereas equation (81) is more useful for oomputatitmal purp<»es. 

Pines (ref. 1) points out that if equation (Bl) is written as 


P__ (p) = cos 
nm 


m 




(B2) 


where A (u) is defined in terms of the Leffendre polynomials P (u) obtained 
nm n 

when m = 0 in equation (Bl) , the following e3q>ressicm results. 







Then , by writing 





(B3) 


(B4) 


equation (96) may be written as 


V = 


00 


A 


n 

p„ A„„(p) cos(m>.)cos^q) + sin(mX)cos™ 9 | 
•^n nm nm nm 

m=0 


(B5) 
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From equation (BS) . Pines recognizes that the terms coe(mX)coe"'q> and 

sin(mX)cos"^(p are merely the real and imaginary parts of the complex number 

(s it)’"'* , where s and t are the direction cosines (x/r) and (y/r ) , respective- 
ly , and i = Equation (96) may be written as 

°° n 

n=0 m=0 

where the mass coefficient function 


= C R„(8.t) + S I„(s,t) 
nm nm m nm m 


(B7) 


and R^(s,t) and are the real and imaginary parts of the complex vari- 

able (s-t-it)^ mentioned previously. That s = cosXoos<p and t = sin X cos q> 
may be established from the right spherical triangle that results from the definition 
of <p and X. 

Equaticm (B6) is an expression for the gravitational potential in terms of the 
magnitude of the position vector and its direction cosines . Its gradient may be 
written as 


vv = 


(i 


av 

5F 


sdV 

r sir 


r at r dpi — 


^ 1 av 


1 av 




1 av 

r 57T 


(B8) 


where i, and k are the unit base vectors along the Cartesian coordinate axes . 
7 is thi gpradient vector operator , and 

R = si + ^ + pk (B9) 

The comparable equation for the gradient of equation (96) is 

= If * sTarvT I f 
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Differentiation of equation (B6) requires differentiation of equaticm (B7) , 
which gives 

dD„„(s.t) dR„(s.t) ai„(8.t) 

nm _ m _ m 

5s '^nm 5s ^nm Tfi 

I nm m-i nm m~i j 

= mE^^(s,t) (BID 

nm 

dD„^(s.t) 3R„(s,t) 

3t nm 5t ®nm ^t 

= m [-C I - (s,t) + S R - (8,t)l 
1^ nm m-i nm m-l j 

= mF^(s,t) (B12) 

since from the Cauchy-Riemann (xxiditions , 

aR„(s,t) ai„(s,t) 
ds 5t 

= mRm_i(s,t) (B13) 

dR„,(8,t) 

= mlm_i(s,t) (B14) 


In equations (Bll) and (B12), respectively, the functions E _(s,t) and F _(s,t) 

nm nm 

are mass coefficient functions. Equations (B13) and (B14) show another advantage 
to the Pines formulation in that these differentiations do not have to be performed 
because the derivatives can be obtained recursively from previous values of the 
real and imaginary parts . This relationship may be seen by writing the complex 
number 


W"™ = Rjjj(s,t) +ilnj(8.t) 


(B15) 
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or 







(B16) 

m ml ml 

(B17) 

I (s.t) = sR ,(s,t) + tR ,(s,t) 
m m-i m~l 

(B18) 


The recursive computation can easily be performed by setting R.fs.t) = 1 and 

lQ(S,t)=0. " 

Differentiation of equation (B6) also requires differentiation of A (p). 
Clearly. 






= A 


h,m+l 


(ii) 


(B19) 


To compute the matrix A^(p) recursively. Pines makes use of the recursive re- 
lations satisfied by all Legendre polynomials. (Spencer (ref. 2) presents five of 
these recursive relations . ) 


^ ['nH = * »'’«<•*> ®20) 

I* ®«> 
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Writing equations (B20) and (B21) in terms of the notation gives 


nAni(|i) - j(ji) = 

Differentiating equations (B22) and (B23) m times gives 


Equation (B25) is the recursive relationship used by Pines in computing the matrix 
Anm(^) • Spencer shows five different relations by which this computation can be 

accomplished. Both authors make use of the initial relationships 

Ann (it) = 1‘3‘5«7‘. , .«(2n - 1) (B26) 

and 

Equation (B27) is obtained from equation (B19) . Equation (B22) is useful in col- 
lecting terms in equation (B8) . Equations (B25) . (B26) . and (B27) reveal another 
desirable feature of this formulation in that errors in the direction cosine p are 
diminished by l/(n - m) . Equation (B8) may now be written as 

VV = + sa^^ i + ^ag + ta^^i + ^ag + pa^^k (B28) 
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where the coefficients are defined as 


n=0 ^ m=0 

oo n 

•2 = E ^ E 

n=0 m=0 


“3" 2 

n=0 * m=0 

•• n 

•4' - E ^ E 

n=0 ^ m=0 

so that 


[®nm^®’*^i 

n=0 m=0 * 

-Vl.m*l<>‘>“nm<»''>£j 

where , from equation (B4) . 

^^n _ ^n+1 


(B29) 


(B30) 


(B31) 


(B32) 


(B33) 


(B34) 
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Equation (B33) is suitable for obtaining higher derivatives of the potential 
function to be used in estunaiing errors realized fircm variations in the harmonic 
coefficients. Pines gives these relations in reference 1, but Spencer (ref. 2) shows 
this work in better form in that he corrects a few errors found in this part of ref- 
erence 1 . Spencer's p^>er is by far the roost complete and accurate , as well as the 
most comprehensible . and should be preferred to reference 1 . 

Sample computer program listings for the computation of the gravitational 
potential functiem and the components of its gradient are provided in the following 
pages of this appendix . The listings are presented primarily to emphasize the 
brevity . 
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